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ABSTRACT 

We present a new algorithm for isolating the real roots of 
a system of multivariate polynomials, given in the mono- 
mial basis. It is inspired by existing subdivision methods 
in the Bernstein basis; it can be seen as generalization of 
the univariate continued fraction algorithm or alternatively 
as a fully analog of Bernstein subdivision in the monomial 
basis. The representation of the subdivided domains is done 
through homographies, which allows us to use only integer 
arithmetic and to treat efficiently unbounded regions. We 
use univariate bounding functions, projection and precon- 
ditionning techniques to reduce the domain of search. The 
resulting boxes have optimized rational coordinates, corre- 
sponding to the first terms of the continued fraction expan- 
sion of the real roots. An extension of Vincent's theorem to 
multivariate polynomials is proved and used for the termina- 
tion of the algorithm. New complexity bounds are provided 
for a simplified version of the algorithm. Examples com- 
puted with a preliminary CH — h implementation illustrate 
the approach. 

Categories and Subject Descriptors 

1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulation — algebraic algorithms 

General Terms 

Algorithms, Theory 

Keywords 

subdivision algorithm, homography, tensor monomial basis, 
continued fractions, C++ implementation 

1. INTRODUCTION 

The problem of computing roots of univariate polynomi- 
als has a long mathematical history [14]. Recently, some 
new investigations focused on subdivision methods, where 
root localization is based on simple tests such as Descartes' 
Rule of Signs and its variant in the Bernstein basis [131 [7J 
[4]. Complexity analysis was developed for univariate integer 
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polynomial taking into account the bitsize of the coefficients, 
and providing a good understanding of their behavior from a 
theoretical and practical point of view. Approximation and 
bounding techniques have been developed [2] to improve the 
local speed of convergence to the roots. 

Even more recently a new attention has been given to 
continued fraction algorithms (CF), see e.g. [161 118] and 
references therein. They differ from previous subdivision- 
based algorithms in that instead of bisecting a given ini- 
tial interval and thus producing a binary expansion of the 
real roots, they compute continued fraction expansions of 
these roots. The algorithm relies heavily on computations 
of lower bounds of the positive real roots, and different ways 
of computing such bounds lead to different variants of the 
algorithm. The best known worst-case complexity of CF 
is Ob(o1 5 t 2 ) [16] . while its average complexity is Ob{o1 3 t), 
thus being the only complexity result that matches, even in 
the average the complexity bounds of numerical algorithms 
[15] . Moreover, the algorithm seems to be the most efficient 
in practice [51 118]. 

Subdivision methods for the approximation of isolated 
roots of multivariate systems are also investigated but their 
analysis is much less advanced. In [17] . the authors used 
tensor product representation in Bernstein basis and domain 
reduction techniques based on the convex hull property to 
speed up the convergence and reduce the number of sub- 
divisions. In [5], the emphasis is put on the subdivision 
process, and stopping criterion based on the normal cone to 
the surface patch. In [12], this approach has been improved 
by introducing pre-conditioning and univariate-solver steps. 
The complexity of the method is also analyzed in terms of 
intrinsic differential invariants. 

This work is in the spirit of [12]. The novelty of our ap- 
proach is the presentation of a tensor-monomial basis al- 
gorithm that generalizes the univariate continued fraction 
algorithm and does not assume generic position. We apply 
a subdivision approach also exploiting certain properties of 
the Bernstein polynomial representation, even though no 
basis conversion takes place. 

Our contributions are as follows. We propose a new adap- 
tive algorithm for polynomial system real solving that acts 
in monomial basis, and exploits the continued fraction ex- 
pansion of (the coordinates of) the real roots. This yields the 
best rational approximation of the real roots. All computa- 
tions are performed with integers, thus this is a division-free 
algorithm. We propose a first step towards the generaliza- 
tion of Vincent's theorem to the multivariate case (Th. I4.2( l 
We perform a (bit) complexity analysis of the algorithm, 



when oracles for lower bounds and counting the real roots 
are available (Prop.[5]2} and we propose non-trivial improve- 
ments for reducing the total complexity even more (Sec, 15. 3)) , 
In all cases the bounds that we derive for the multivariate 
case, match the best known ones for the univariate case, if 
we restrict ourselves to n = 1. 

1.1 Notation 

For a polynomial / g M{xx, .., x n ], deg(/) denotes its to- 
tal degree, while deg x .(f) denotes its degree w.r.t. Xi. Let 
/fe) = f(xi,..,x n ) G R[xi,..,x n ] with deg Xk f = d k , k = 
l,..,n. If not specified, we denote d = d(f) — max{di, .., d n }. 

We are interested in isolating the real roots of a system 
of polynomials fi(x), .., f s (x) G Z[xx, .., x n ], in a box 7o = 
[ui,«i] x •■• x [un,v n ] C R", u k ,v k G Q. We denote by 
•2k»(/) ={?6 K n ; f(p) = 0} the solution set in K" of the 
equation f(x) = 0, where K is R or C. 

In what follows Ob, resp. O, means bit, resp. arithmetic, 
complexity and the Ob, resp. O, notation means that we 
are ignoring logarithmic factors. For a G Q, C (a) > 1 is the 
maximum bit size of the numerator and the denominator. 
For a polynomial / g Z[xx, .., x n ], we denote by C (/) the 
maximum of the bitsize of its coefficients (including a bit 
for the sign). In the following, we will consider classes of 
polynomials such that log(d(/)) = 0(C (/))• 

Also, to simplify the notation we introduce multi-indices, 
for the variable vector x — (xx, ..,x n ), x L := x 1 ^ ■ ■ ■ x 1 ™ , the 

1 -eT- (f) - (t)- (t)- - 

tensor Bernstein basis polynomials of multidegree degree d 
of a box / are denoted B(x;i,d; I) := (xx', Ux, Ux) • • ■ 
B l d " n (x n ;u n ,Un) where / = [u,v] := [tti,«i] x ■■• x [u n ,v n ]. 

1.2 The general scheme 

In this section, we describe the family of algorithms that 
we consider. The main ingredients are 

• a suitable representation of the equations in a given 
(usually rectangular) domain, for instance a represen- 
tation in the Bernstein basis or in the monomial basis; 

• an algorithm to split the representation into smaller 
sub- domains; 

• a reduction procedure to shrink the domain. 

Different choices for each of these ingredients lead to algo- 
rithms with different practical behaviors. The general pro- 
cess is summarized in Alg. [T] 

The instance of this general scheme that we obtain gen- 
eralizes the continued fraction method for univariate poly- 
nomials; the realization of the main steps (b-e) can be sum- 
marized as follows: 

b) Perform a precondition process and compute a lower 
bound on the roots of the system, in order to reduce 
the domain. 

c) Apply interval analysis or sign inspection to identify 
if some fi has constant sign in the domain, i.e. if the 
domain contains no roots. 

d) Apply Miranda test to identify if the domain contains 
a single root. In this case output (I, fx, .., fs)- 



Algorithm 1.1: Subdivision scheme 



Input: A set of equations fx, fz, ■■, fs G Z[^] 

represented over a domain /. 
Output: A list of disjoint domains, each containing one 

and only one real root of fx = • • • = f s = 0. 
Initialize a stack Q and add (/, fx, .., f s ) on top of it; 
While Q is not empty do 

a) Pop a system (7, fx, .., f s ) and: 

b) Perform a precondition process and/or a reduction 
process to refine the domain. 

c) Apply an exclusion test to identify if the domain 
contains no roots. 

d) Apply an inclusion test to identify if the domain 
contains a single root. In this case output 
(I,fi, ..,/.). 

e) If both tests fail split the representation into a 
number of sub-domains and push them to Q. 



e) If both tests fail, split the representation at (1,..,1) 
and continue. 

In the following sections, we are going to describe more 
precisely the specific steps and analyze their complexity. In 
Sec. [2] we describe the representation of domains via ho- 
mographies and the connection with the Bernstein basis 
representation. Subdivision, based on shifts of univariate 
polynomials, reduction and preconditionning are analyzed 
in Sec. [3] Exclusion and inclusion tests as well as a gener- 
alization of Vincent's theorem to multivariate polynomials, 
are presented in Sec. [4] In Sec. we recall the main proper- 
ties of Continued Fraction expansion of real numbers and use 
them to analyze the complexity of a subdivision algorithm 
following this generic scheme. We conclude with examples 
produced by our C++ implementation in Sect. [6] 

2. REPRESENTATION: HOMOGRAPHIES 

A widely used representation of a polynomial / in a box 
is the tensor-Bernstein representation. De Casteljau's algo- 
rithm provides an efficient way to split this representation to 
smaller domains. A disadvantage is that converting integer 
polynomials to Bernstein form results in rational Bernstein 
coefficients. We follow an alternative approach that does 
not require basis conversion since it applies to monomial 
basis: We introduce a tensor-monomial representation, i.e. 
a representation in the monomial basis over P 1 x ■ ■ ■ x P 1 
and provide an algorithm to subdivide this representation 
analogously to the Bernstein case. 

In a tensor-monomial representation a polynomial is rep- 
resented as a tensor (higher dimensional matrix) of coeffi- 
cients in the natural monomial basis, that is, 

di,..,d n d 

for every equation / of the system. Splitting this represen- 
tation is done using homographies. The main operation in 
this computation is the Taylor shift. 



Definition 2.1. A homography (or Mobius transforma- 
tion) is a bijective projective transformation Tri = (Hi, .., Tin) 
defined over P 1 x ■ ■ • x P 1 as 



x k i-> H k (x k ) 



Ct k X k + Pk 



^ k x k + S k 

with Qfc,/3fc,7fe,<5 fe G Z, 7fc(5 fe t^O, k = 1, „,n. 

Using simple calculations, we can see that the inverse 

8 k x k - (3 k 



w* 0=*) 



is also a homography. Also, notice that if detTi > then, 
taking proper limits when needed, we can write 



n 
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hence H(f) 



H(f) ~\{{lkXk+5k) dk -{fon){x) 



is a polynomial defined over R™ that corresponds to the 
(possibly unbounded) box 
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of the initial system, in the sense that the zeros of the ini- 
tial system in Ih are in one-to-one correspondence with the 
positive zeros of H(f). 

We focus on the computation of H(f). We use the base 
homographies T k (f) = f\ Xk =x k +c (translation by c) or sim- 
ply T k (f) if c = 1, Cfc(/) = f\x h =cx k (contraction by c) and 
Rk(f) = Xu k f\ Xh= i/ x . (reciprocal polynomial) . These nota- 
tions are naturally extended to variable vectors; for instance 
T- = (T 1 G1 ,..,T 1 J n ), c = (ci,..,Cn) 6 Z n . Complexity results 
for these computations appear in the following sections. We 
can see that they suffice to compute any homography: 

Lemma 2.2. The group of homographies H with coeffi- 
cients in Z is generated by Rk,C k ,T k , k = 1, ..,n, c G Z. 

Proof. It can be verified that Hk(f) is computed as 



HkU) = Cl k RkC & k k TkRkC l 



l 3 k/ S k- a k/lkrpl 3 k/ S k 



'(/) 



where the product denotes composition. We abbreviate 



l/c 



RkC k Rk and T, 



l/c 



C%TlR k T£R k . □ 



Representation via homography is not far from Bernstein 
representation: 

Lemma 2.3. Let f = X^Lo ^ B"(x, Ih) the Bernstein ex- 
pansion of f in the box Ih yielded by a homography H. If 

d 
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Proof. Let [u k ,v k ] = [^fif^]- For a tensor-Bernstein 

polynomial [ ~ ] - — - — r-r(x — uY(v — x) d ~ l we compute 
\i I (v — u) d 
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as needed. □ 

Corollary 2.4. The Bernstein expansion of f in Ih is 



E 



-Bi (x;I H ). 



That is, the coefficients of H(f) coincide with the Bernstein 
coefficients up to contraction and binomial factors. 

Thus tensor-Bernstein coefficients and tensor-monomial 
coefficients in a sub-domain of differ only by multipli- 
cation by positive constant. In particular they are of the 
same sign. Hence this corollary allows us to take advantage 
of sign properties (eg. the variation diminishing property) 
of the Bernstein basis without computing it. 

The resulting representation of the system consists of the 
transformed polynomials H(fi), ..,H(f n ), represented as ten- 
sors of coefficients as well as 4n integers, a k , /3fe,7fe, Sk for 
k — 1, .., n from which we can recover the endpoints of the 
domain, using ([2]). 

3. SUBDIVISION AND REDUCTION 
3.1 The subdivision step 

We describe the subdivision step using the homography 
representation. This is done at a point u — (ui,..,u„) G 
Z> . It consists in computing up to 2 n new sub-domains 
(depending on the number of nonzero W^'s), each one having 
u as a vertex. 

Given H(fi), .., H(f s ) that represent the initial system at 
some domain, we consider the partition of R™ defined by 
the hyperplanes x k = u k , k — l,..,n. These intersect at 
u hence we call this partition at u. Subdividing at u is 
equivalent to subdividing the initial domain into boxes that 
share the common vertex Tt(u) and have faces either parallel 
or perpendicular to those of the initial domain. 

We need to compute a homography representation for ev- 
ery domain in this partition. The computation is done co- 
ordinate wise; observe that for any domain in this partition 
we have, for all k, either Xk G [0, u k ] or Xk G [life, 00]. It 
suffices to apply a transformation that takes these domains 
to R+. In the former case, we apply TlR k C^ k to the current 
polynomials and in the latter case we shift them by u k , i.e. 
we apply T^ k . The integers a k , f3 k ,j k , 8 k that keep track of 
the current domain can be easily updated to correspond to 
the new subdomain. 



We can make this process explicit in general dimension: 
every computed subdomain corresponds to a binary number 
of length n, where the k— th bit is 1 if T k RkC Uk is applied 
or if T^ k is applied. 

In our continued fraction algorithm the subdivision is per- 
formed at u = 1. 

Illustration. Let us illustrate this process in dimension 
two. The system is defined over R>o. We subdi- 

vide this domain into [0, l] 2 , [0, 1] x R>i, K>i x [0, 1] and 
R>i x R>i . Equivalently, we compute four new pairs of poly- 
nomials, as illustrated in Fig.[T](we abbreviate S k = Tj~R k ). 



(0,0) 



SiT 2 {f) 


Z\T 2 (/) 




(1,1) 











Figure 1: Subdividing the domain of /. 

Complexity of subdivision step. The transformation of 
a polynomial into two sub-domains, i.e. splitting w.r.t. one 
direction, consists of performing d"" 1 univariate shifts, one 
for every coefficient 6 Z[xk] of / £ Z[a;fc][a;i, ..,Xk, ..,!„]. 

If the subdivision is performed in every direction, each 
transformation consists of rf n_1 univariate shifts for every 
variable, i.e. nd n ~ shifts. There are 2™ sub-domains to 
compute, hence a total of n 2 2 n d n ~ 1 shifts have to be per- 
formed in a single subdivision step. We must also take into 
account that every time a univariate shift is performed, the 
coefficient bitsize increases. 

The operations 

Tk(f) = /k=* fc +i and T h R k (f) = (x k + l) d *f\ i 

are essentially of the same complexity, except that the sec- 
ond requires one to exchange the coefficient of c^,..,^,..,^ 
with Ci x$ ..,d k -i k ,..,i n before translation, i.e. an additional 
0(d n ) cost. Hence we only need to consider the case of 
shifts for the complexity. 

The continued fraction algorithm subdivides a domain us- 
ing unit shifts and inversion. Successive operations of this 
kind increase the bitsize equivalently to a big shift by the 
sum of these units. Thus it suffices to consider the general 
computation of f(x + u) to estimate the complexity of the 
subdivision step. 

Lemma 3.1 (Shift complexity). The computation of 
f{x_ + u) with £.(/) = t and C{u k ) < a, k = 1, .., n can be 
performed in OB(n 2 d n r + d n+1 n 3 a). 

Proof. We use known facts for the computation of T^ h (/) 
for univariate polynomials. If deg fc / = dk and / is univari- 
ate, this operation is performed in OB(d k a+dkT); the result- 
ing coefficients are of bitsize T+dk<7 [20]. Hence f{x\, .., Xk + 
u k ,..,x n ) is computed in Os(d' I_1 (d^cr + d k T)). 

Suppose we have computed f(xi+ui, x k -i+u k -i, x k , 
for some k. The coefficients are of bitsize r + "Y^Zl o%- The 
computation of shift w.r.t. k— th variable f(x\ +ui, ..,Xk + 



Uh,Xk+i, ■■,x n ) results in a polynomial of bitsize t+X).* o~i 
and consists of d n_1 Os(d 2 X^=i °« + ^ r )) operations. That 
is, we perform d n_1 univariate polynomial shifts, one for 
every coefficient of / in Z[xk]\xi, ..,Xk, .., x n ]. 

This gives a total cost for computing f(x_ + u) of 

n / k \ n 

d n_1 ( d 2 en + dr J = nd n r + d n+1 ^^(w + 1 — k)ok- 

fe=l V i=l / fc=l 

The latter sum implies that it is faster to apply the shifts 
with increasing order, starting with the smallest number Uk- 
Since dk = O(o) for all k, and we must shift a system of 
0(n) polynomials we obtain the stated result. □ 

Let us present an alternative way to compute a sub-domain 
using contraction, preferable when the bitsize of u is big. 
The idea behind this is the fact that T£ and T^C C compute 
the same sub-domain, in two different ways. 

Lemma 3.2. /// = J2f= c i2r< c (f) = T > then the coeffi- 
cients ofC u (f), C(uk) < o~, k — 1, ..,n, can be computed in 
B (d n r + nd n+1 a) . 

Proof. The operation, i.e. computing the new coeffi- 
cients Civr- can be done with 0(d n ) multiplications: Since 
y(ii,..,i*.—»n) _ Ufc ^(*i,..,ih-i,..,»n) j jf these powers are com- 
puted successively then every coefficient is computed using 
two multiplications. Moreover, it suffices to keep in memory 
the n powers u ( t i,-,»f C -i> i fe- 1 > ! fe+i,,-,»n) j = i } _, >n i n order 
to compute any M-Cj. Geometrically this can be understood 
as a stencil of n points that sweeps the coefficient tensor 
and updates every element using one neighbor at each time. 
The bitsize of the multiplied numbers is 0(t + da) hence 
the result follows. □ 

Now if we consider a contraction followed by a shift by 
1 w.r.t. Xk for 0(n) polynomials we obtain Os(n 2 d n r + 
n 3 d n+1 + nd n+1 a) operations for the computation of the 
domain. The disadvantage is that the resulting coefficients 
are of bitsize 0(t + da) instead of 0(t + na) with the use 
of shifts. Also note that this operation would compute a ex- 
pansion of the real root which differs from continued fraction 
expansion. 

3.2 Reduction: Bounds on the range of / 

In this section we define univariate polynomials whose 
graph in R n+1 bounds the graph of /. For every direction 
k, we provide two polynomials bounding the values of / in 
K" from below and above respectively. 

Define 

m k (f;x k ) = ^2 "ii 11 c n--i n x l k ( 3 ) 
d k 

Mk(f;x k ) = 22 max c il .. in x l k k (4) 

Lemma 3.3. For any x £ R", n > 1 and any k = 1, ..,n, 
we have 

m k {f; x k ) < ^ < M k (f; x k ). (5) 

s^k i B =0 



Proof. For x G R™ , we can directly write 

(d k \ d s 

max c n ..i„ a# ) II ]C x * a ' 
i h =0 n ' / s#fci s =0 

The product of power sums is greater than 1; divide both 
sides by it. Analogously for M k (f; x k ). □ 

Corollary 3.4. Given k G {1, .., n}, if u & R™ itfiin 
u fc G]0, /Life], w/iere 

fmin. pos. roo£ of M k (f,x k ) if Mfc(/;0) < 
nun. pos. root of m k ( f,x k ) i/ra fe (/;0)>0 , 
otherwise 

then f(u) 7^ 0. Consequently, all positive roots of f lie in 
R> M1 x ■ • ■ x R> M „ . .4/so, /or it 6 R+ u fe G LMs, oo], 

{max. pos. root of Mk(f,Xk) if Mk(/;oo) < 
rnaa;. pos. root of m k (f,x k ) ifm k (f;oo) >0 , 
oo otherwise 

it is f{u) 7^ 0, i.e. a// pos. roots are in R<Aii X • • • X R<vw„ ■ 
Combining both bounds we deduce that x ••• x 

\p n ,M„] is a bounding box for / _1 ({0}) n R+. 

Proof. The denominator in ([5]) is always positive in R™ . 
Let w g 1" with u k G [0,/j ft ]. If M k (f,0) < then also 
M k (f,u) < and it follows f(u) < 0. Similarly m fc (/, 0) > 
=> m k (f,u) > => /(m) < 0. The same arguments hold 
for [M k ,oo], M fc (/;oo) = fi(M fc (/;a; fc ))(0), m fc (/;oo) = 
R{ m k{f', Xk))(Q), and R(f), since lower bounds on the zeros 
of yield upper bounds on the zeros of /. □ 

Thus lower and upper bounds on the k— th coordinates of 
the roots of (/i, .., f s ) axe given by 

max {Hk(fi)} and min {M k (fi)} (6) 

i=\, .. ,s i — 1, . . ,s 

respectively, i.e. the intersection of these bounding boxes. 

We would like to remain in the ring of integers all along 
the process, thus integer lower or upper bounds will be 
used. These can be the floor or ceil of the above roots of uni- 
variate polynomials, or even known bounds for them, e.g. 
Cauchy's bound. 

If the minimum and maximum are taken with the order- 
ing of coefficients defined as c; -< Cj -<=>■ <k(~)'f 8 d ~^ < 

Cj( d ^ l 8 d ~ l then different m k (f,x k ),M k (f,x k ) polynomials 
are obtained. By Cor. I2.4l their control polygon is the lower 
and upper hull respectively of the projections of the tensor- 
Bernstein coefficients to the fc— th direction and are known 
to converge quadratically to simple roots when precondi- 
tioning (described in the following paragraph) is utilized |121 
Cor. 5.3]. 

Complexity analysis. The analysis of the subdivision step 
in Sect. 13.21 applies as well to the reduction step, since re- 
ducing the domain means computing a new subdomain and 
ignoring the remaining part. 

If a lower bound / is known, with C(l k ) = O(a), then the 
reduction step is performed in OB{n 2 d n r + d n+1 n 3 a). This 
is an instance of Lem. 13.11 

The projections of Lem. 13.31 are computed using 0(d n ) 
comparisons. The computation of I costs OB(d S T) in aver- 
age, for solving these projections using univariate CF algo- 
rithm. Another option would to compute well known lower 
bounds on their roots, for instance Cauchy's bound in 0(d). 



Illustration. Consider a bi-quadratic /o G R[x,j/], namely, 
deg^/o = deg X2 fo — 2 with coefficients dj. Suppose that 
fo — H(f) for Io = Ih- We compute 

2 2 

miff, xi) = > min ca x 1 and Mi fi) = > max cax 1 . 

^ — ' 7—0, ..2 t— 1 1=0, ..2 

i=0 i=0 

thus mix) < ^ 1, T 2 \ < M(x). Fig. shows how these 
univariate quadratics bound the graph of / in Iq. 




Figure 2: The enveloping polynomials M\[x), m\(x) 
in domain Jo for a bi-quadratic polynomial f(x,y). 

3.3 Preconditioning 

To improve the reduction step, we use preconditioning. 
The aim of a preconditioner is to tune the system so that 
it can be tackled more efficiently; in our case we aim at 
improving the bounds of Cor. 13.41 

A preconditioning matrix P is an invertible s x s matrix 
that transforms a system (/i, .., / s )' into the equivalent one 
P ' (/li - J fs) ■ This transformation does not alter the roots 
of the system, since the computed equations generate the 
same ideal. The bounds obtained on the resulting system 
can be used directly to reduce the domain of the equations 
before preconditioning. Preconditioning can be performed 
to a subset of the equations. 

Since we use a reduction process using Cor. 13.41 we want 
to have among our equations n of them whose zero locus 
/ _1 ({0}) is orthogonal to the k— th direction, for all k. 

Assuming a square system, we precondition H(f\), .., H( f n ) 
to obtain a locally orthogonal to the axis system; an ideal 
preconditioner would be the Jacobian of the system evalu- 
ated at a common root; instead, we evaluate Jn(f) at the im- 
age of the center u of the initial domain Ih, u k = ak5k+ P k ~> k 
Thus we must compute the inverse of the Jacobian matrix 
Ja(f)(x) = [d Xi H(fj)(x)]i<ij } < n evaluated at u := H(u) = 
{Si/ji, ..,<5 n /7„). 

Precondition step complexity. Computing Jh{J){u) ■ 
(H(fi), ..,H(fn)Y is done with cost Os(n 2 d n ) and evaluat- 
ing at u has cost Os(n 2 d n_1 ). We also need Ob(u 2 ) for 
inversion and 0(n 2 d n ) for multiplying polynomials times 
scalar as well as summing polynomials. This gives a precon- 
dition cost of order 0(n 2 d n ). 

4. EXCLUSION - INCLUSION CRITERIA 

A subdivision scheme is able to work when two tests are 
available: one that identifies empty domains (exclusion test) 
and one that identifies domains with exactly one zero (inclu- 
sion test). If these two tests are negative, a domain cannot 



be neither included nor excluded so we need to apply fur- 
ther reduction/subdivision steps to it. The certification is 
the following: if the result of the test is affirmative, then 
this is undoubtedly true. 

Exclusion test. The bounding functions defined in the pre- 
vious section provide a fast filter to exclude empty domains. 
Define min{} = oo and max{} = 0. 

Corollary 4.1. If for some k £ {1, ..,n} and for some 
i 6 {1, .., s} it is Hk(fi) = oo or Mk(fi) = then the system 
has no solutions. Also, if maxi = i i .. iS {//fe(/i)} > minj = i,.. lS 
{Mlk(fi)} then there can be no solution to the system. 

Proof. For the former statement observe that fi has no 
real positive roots, thus the system has no roots. The latter 
statement means that the reduced domains of each fi, i = 
1, .., s do not intersect, thus there are no solutions. □ 

We can use interval arithmetic to identify additional empty 
domains; if the sign of some initial fi is constant in Ih = 
7i(R™ ) then this domain is discarded. We can also simply 
inspect the coefficients of each H(fi); if there are no sign 
changes then there corresponding box contains no solution. 

The accuracy of these criteria greatly affects the perfor- 
mance of the algorithm. In particular, the sooner an empty 
domain is rejected the less subdivisions will take place and 
the process will terminate faster. We justify that the exclu- 
sion criteria will eventually succeed on an empty domain by 
proving a generalization of Vincent's theorem to the tensor 
multivariate case. 

Theorem 4.2. Let f(x) = J2f=o C *-— L be a polynomial 
with real coefficients, such that it has no (complex) solu- 
tion with ^t(zk) > for k = 1, .., n. Then all its coefficients 
Ci 1 ,..,i n are of the same sign. 

Proof. We prove the result by induction on n, the num- 
ber of variables. For n — 1, this is the classical Vincent's 
theorem [T]. 

Consider now a polynomial 

f(xi,x 2 ) = ^2 Ci ui2 xl 1 x\ 2 

0<i 1 <d ± ,0<i 2 <d 2 

in two variables with no (complex) solution such that ^R(xi) > 
for i — 1,2. We prove the result for n = 2, by induction 
on the degree d = d 1 + 0I2 . The property is obvious for poly- 
nomials of degree d = 0. Let us assume it for polynomials 
of degree less than d. 

By hypothesis, for any Zj £ C with > 0, the uni- 

variate polynomial f(z\,X2) has no root with 5R(x2) > 0. Ac- 
cording to Lucas theorem the complex roots of d X2 f(zi, X2) 
are in the convex hull of the complex roots of f(zi,X2). 
Thus, there is no root of d X2 f(x\,X2) with ^R{x\) > and 
$l(x2) > 0. By induction hypothesis, the coefficients of 
d X2 f(xi,X2) are of the same sign. We decompose P as 

f(x!,x 2 ) = f(xi,0) + f 1 (x 1 ,x 2 ) 

where fi(x 1 ,x 2 ) = Eckh^!,!^^, c *i>*2 x i x 2 with c H , i2 
of the same sign, say positive. By Vincent theorem in one 
variable, as /(a;i,0) has no root with 5R(a:i) > 0, the coeffi- 
cients 0^,0 of /(a;i,0) are also of the same sign. If this sign 
is different from the sign of Ci lt i 2 for 12 1 (ie. negative 
here), then /(0, X2) has one sign variation in its coefficients 
list. By Descartes rule, it has one real positive root, which 



contradicts the hypothesis on /. Thus, all the coefficients 
have the same sign. 

Assume that the property has been proved for polyno- 
mials in n — 1 variables and let us consider a polynomial 
fis) = 53*=o Ci — 1X1 n var i aD l es with no (complex) solution 
such that $l(xk) > for k — l,..,n. For any zi,..,z n -i € 
C with 5R(zfc) > 0, for k = l,..,n — 1, the polynomial 
f(zi,.., z n -\,x n ) and d x „f(zi,.., z n -i,x n ) has no root with 
5ft(x n ) > 0. By Lucas theorem and induction hypothesis on 
the degree, d Xn f(x) has coefficients of the same sign. We 
also have f(x\, ..,x n -i, 0) with coefficients of the same sign, 
by induction hypothesis on the number of variables. If the 
two signs are different, then /(0, .., 0, x n ) has one sign vari- 
ation in its coefficients and thus one real positive root, say 
Ca, which cannot be the case, since (0, ..,0,£ n ) would yield 
a real root of /. We deduce that all the coefficients of / are 
of the same sign. 

This completes the induction proof of the theorem. □ 

This implies that empty regions will be eventually ex- 
cluded by sign inspection. 

Corollary 4.3. Let H(f) = J2f= c lS. l be the 

represen- 
tation of f through Ti in a box Ih = \u, v\ ■ If there is no 
toot z £ C" of f such that 



then all the coefficients a 1 ,..,i n are of the same sign. 

That is, if dist 00 (Zc"(f),m) > e, where m is the center 
of Ih, then Ih is excluded by sign conditions. 

Proof. The interval [uk,Vk] is transformed by HT 1 into 
[0, +00] and the disk \z h - u *+"* | < " fc ~" fc is transformed 
into the half complex plane 5ft(zfc) > 0. We deduce that 
H(f) has no root with $l(z k ) >0,k = l, ..,n. By Thm.l4~2l 
the coefficients of H(f) are of the same sign. □ 

We deduce that if a domain is far enough from the zero 
locus of some fi then it will be excluded, hence redundant 
empty domains concentrate only in a neighborhood of / = 0. 

Definition 4.4. The tubular neighborhood of size e of fi 
is the set 

r E (/,) = (i6l" : 3z£C",/,(z)=0, si \\z - x\\^ < e}. 

We bound the number of boxes that are not excluded at 
each level of the subdivision tree. 

Lemma 4.5. Assume that for so > 0, C\iT Eo (fi) n Io is 
bounded. Then the number of boxes of size e < £0 kept by 
the algorithm is less than (1 + -^) n c, where c > is such 
that Ve st. £0 > £ > 0, 

V(f, e) := volume (n| =1 r E (/ 4 ) n Io)) < ce n . 

Proof. Consider a subdivision of a domain Io into boxes 
of size £ < eq. We will bound the number N of boxes in 
this subdivision that are not rejected by the algorithm. By 
Cor. 14.31 if a box is not rejected, then we have for all i — 
1, ..,s distoo(2c»(/i), m) < £, where m is the center of the 
box. Thus all the points of this box are at distance < e(l + 
^) to 2fc»(/<) that is in nf =1 r (fi) n I . 

To bound N, it suffices to estimate the n— dimensional 
volume V(f, e), since we have: 

Ne n < volume (^^j (/*) n Jo) =V(f,e(l+^)). 



When e tends to 0, this volume becomes equivalent to a 
constant times e n . For a square system with single roots 
in Jo, it becomes equivalent to the sum for all real roots 
£ in Io of the volumes of parallelotopes in n dimensions 
of height 2e and unitary edges proportional to the gradi- 
ents of the polynomials evaluated at the common root; It is 

thus bounded by e n 2 n X}<e/ fTTTvTT^rn ' ^ e deduce that 

there exists a constant c > 2 n Y]*^ , ,-, i'f/ffl,L M such that 
— ^Ceio Hi l|v/ 4 (C)|| 

V(f, e) < ce n < oo. For overdetermined systems, the vol- 
ume is bounded by a similar expression. Since V(f, e)e~" 
has a limit when e tends to 0, we deduce the existence of the 
finite constant c and the bound of the lemma on the number 
of kept boxes of size e. □ 

Inclusion test. We present a test that discovers common 
solutions, in a box, or equivalently in R™, through homogra- 
phy. To simplify the statements we assume that the system 
is square, i.e. s — n. 

Definition 4.6. The lower face polynomial of f w.r.t. 
direction k is low(/, fc) = f\ Xk=0 . The upper face polynomial 
of f w.r.t. k is upp(/,fe) = ,f\ Xk =oo ■= Rk(.f)\x k =o- 

Lemma 4.7 (Miranda Theorem [21]). If for some per- 
mutation n : {l,..,n} — > {l,..,n}, sign(low(.ff (/*,), 7r(fc))) 
and sign(upp(//(/fe), n(k))) are constant and opposite for all 
k — l,..,n, then the equations (/l, •-,/») have at least one 
root in Ih- 

The implementation of the Miranda test can be done effi- 
ciently if we compute a — 1 matrix with th entry 1 
iff sign (low (//(/;), j)) and sign(upp(i7(/i), j)) are opposite. 
Then, Miranda test is satisfied iff there is no zero row and 
no zero column. To see this observe that the matrix is the 
sum of a permutation matrix and a — 1 matrix iff this 
permutation satisfies Miranda's test. 

Combined with the following simple fact, we have a test 
that identifies boxes with a single root. 

Lemma 4.8. 7/ det Jf(x) has constant sign in a box I, 
then there is at most one root of f = (fi, ..,/n) wi 7. 

Proof. Suppose u,v G I are two distinct roots; by the 
mean value theorem there is a point w on the line segment 
uv, and thus in I, s.t. Jf(w) ■ (u — v) = /(it) — f(v) — 
hence det Jf (w) = 0. □ 

Complexity of the inclusion criteria. Miranda test can 
be decided with 0(n 2 ) evaluations on interval (cf. [§]) as 
well as one evaluation of Jf, overall 0{n 2 d n ) operations. 
The cost of the inclusion test is dominated by the cost of 
evaluating 0(n) polynomials of size 0(d n ) on an interval, 
i.e. 0(nd n ) operations suffice. 

Proposition 4.9. If the real roots of the square system 
in the initial domain Iq are simple, then Alg. \7\ stops with 
boxes isolating the real roots in Io- 

Proof. If the real roots of / = (/i, .., f n ) in Io are sim- 
ple, in a small neighborhood of them the Jacobian of / has 
a constant sign. By the inclusion test, any box included in 
this neighborhood will be output if and only if it contains 
a single root and has no real roots of the jacobian. Other- 
wise, it will be further subdivided or rejected. Suppose that 
the subdivision algorithm does not terminate. Then the size 



of the boxes kept at each step tends to zero. By Cor. 14.31 
these boxes are in the intersection of the tubular neighbor- 
hoods (nf =1 tub e (/i)) n R" for s > the maximal size of 
the kept boxes. If e is small enough, these boxes are in a 
neighborhood of a root in which the Jacobian has a constant 
size, hence the inclusion test will succeed. By the exclusion 
criteria, a box domain is not subdivided indefinitely, but is 
eventually rejected when the coefficients become positive. 
Thus the algorithm either outputs isolating boxes that con- 
tains a real root of the system or rejects empty boxes. This 
shows, by contradiction, the termination of the subdivision 
algorithm. □ 

5. THE COMPLEXITY OF MCF 

In this section we compute a bound on the complexity of 
the algorithm that exploits the continued fraction expansion 
of the real roots of the system. Hereafter, we call this algo- 
rithm MCF (Multivariate Continued Fractions). Since the 
analysis of the reduction steps of Sec. [3] and the Exclusion- 
Inclusion test of Sec. [4] would require much more develop- 
ments, we simplify the situation and analyze a variant of 
this algorithm. We assume that two oracles are available. 
One that computes, exactly, the partial quotients of the pos- 
itive real roots of the system, and one that counts exactly 
the number of real roots of the system inside a hypercube 
in the open positive orthant, namely . In what follows, 
we will assume the cost of the first oracle is bounded by Ci , 
while the cost of the second is bounded by C2, and we derive 
the total complexity of the algorithm with respect to these 
parameters. In any case the number of reduction or subdi- 
vision steps that we derive is a lower bound on the number 
of steps that every variant of the algorithm will perform. 
The next section presents some preliminaries on continued 
fractions, and then we detail the complexity analysis. 

5.1 About continued fractions 

Our presentation follows closely [7] . For additional details 
we refer the reader to, e.g., |22l [3l [19]. In general a simple 
(regular) continued fraction is a (possibly infinite) expression 
of the form 

Co H — = [c ,ci,c 2 , ..], 



where the numbers a are called partial quotients, d £ Z and 
Ci > 1 for i > 0. Notice that Co may have any sign, however, 
in our real root isolation algorithm Co > 0, without loss of 
generality. By considering the recurrent relations 

P-l = 1, Po = CO, Pn+1 = Cn + 1 Pn + Pn-1, 
Q-l = 0, Qo = 1, Qn+l = C n +i Q n + Q„-l, 

it can be shown by induction that R n = ^ = [co, cx, .., c n ], 
for n = 0, 1,2, ... 

If 7 = [c ,ci,..] then 7 = co + ^ - ^ + .. = 

co + S^Lj q~ 1 \q an d since this is a series of decreasing al- 
ternating terms it converges to some real number 7. A finite 
section R n — = [co,Ci, ..,c n ] is called the n— th conver- 
gent (or approximant) of 7 and the tails 7,1+1 = [c n +i, Cn+2, ■■] 
are known as its complete quotients. That is 7 = [co, ci, .., c n , 
Tn+i] for n = 0, 1,2, ... There is an one to one correspon- 
dence between the real numbers and the continued fractions, 



where evidently the finite continued fractions correspond to 
rational numbers. 

It is known that Q n > F n+l and that F n+l < <f> n < 
Fn+2, where F n is the n— th Fibonacci number and </> — 
is the golden ratio. Continued fractions are the best 
rational approximation(for a given denominator size). This 
is as follows: 



Qn(Qn + l + Qr< 



< 



7- 



Pn 

Qn 



< 



<0 



-2n+l 



(7) 



Let 7 = [co,ci,..] be the continued fraction expansion of 
a real number. The Gauss-Kuzmin distribution [3 > states 
that for almost all real numbers 7 (meaning that the set of 
exceptions has Lebesgue measure zero) the probability for a 
positive integer 8 to appear as an element c; in the continued 
fraction expansion of 7 is 



Prob[ci — 8] ^2 lg 



(<5+l) 2 



for any fixed i > 0. 



(8) 



8(8 + 2) ' 

The Gauss-Kuzmin law induces that we can not bound the 
mean value of the partial quotients or in other words that the 
expected value (arithmetic mean) of the partial quotients is 
diverging, i.e. 



E[a] = Y2 S Prob i c i = s ] 



for i > 0. 



Surprisingly enough the geometric (and the harmonic) 
mean is not only asymptotically bounded, but is bounded 
by a constant, for almost all 7 £ R. For the geometric mean 
this is the famous Khintchine's constant [10], i.e. 



lim 



K, = 2.685452001. 



It is not known if K. is a transcendental number. The ex- 
pected value of the bit size of the partial quotients is a con- 
stant for almost all real numbers, when n — > 00 or n suf- 
ficiently big [ID]. Notice that in flHJ, i > 0, thus 7 G E is 
uniformly distributed in (0, 1). Let C (a) = bi, then 



E[bi] = OQgK.) = 0(l). 



O) 



5.2 Complexity results 

Let a be an upper bound on the bitsize of the partial 
quotient that appear during the execution of the algorithm. 

Lemma 5.1. The number of reduction and subdivision steps 
that the algorithm performs is Oijr^Td 2 "^ 1 ) . 

PROOF. Let £ = (Cij-iC™) be a real root of the system. 
It suffices to consider the number of steps needed to isolate 
the i coordinate of £. 

Recall, that we assume, working in the positive orthant, 
we can compute exactly the next partial quotient in each 
coordinate; in other words a vector I = (h,..,l n ), where 
each k, 1 < i < n, is the partial quotient of a coordinate of 
a positive reaQ solution of the system. 

Let fe(£) be the number of steps needed to isolate the i 
coordinate of the real root The analysis is similar to the 
univariate case. The successive approximations of Q by the 

1 Actually the analysis holds even in the case where each U 
is the partial quotient of the positive imaginary part of a 
coordinate of a solution of the system. 



lower bound k, yield the fci(£)-th approximant, 
which using (0 satisfies 

P ki(0 



ofO, 



'MO 



< 



< 



-2MO+1 



In order to isolate it suffices to have 



Q 



MO 



< A,(C), 



where Aj(C) is the local separation bound of that is the 
smallest distance between and all the other i-coordinates 
of the positive real solutions of the system. 

Combining the last two equations, we deduce that to achieve 
the desired approximation, we should have 0- 2 MO+i < 
Ai(C), or fe(C) > I - |lgAi(C). That is to isolate the i 
coordinate it suffices to perform 0{— \ lg Ai(C)) steps. To 
compute the total number of steps, we need to sum over all 
positive real roots and multiply by n, which is the number 
of coordinates, that is 

n HO < n-R-n- lg Ai(C) = n-R~n\ lg ]J A 4 (0, 



CSV 



CSV 



where \ V\ — R is the number of positive real roots. 

To bound the logarithm of the product, we use DMM n [8], 
i.e. aggregate separation bounds for multivariate, zero-di- 
mensional polynomial systems. It holds 



2nrd 2n - 1 + 2nd n lg(nd 2 



Taking into account that R < d n we conclude that the num- 



ber of steps is 0(r?Td 2n -^) 



□ 



Proposition 5.2. The total complexity of the algorithm 
is d B {2 n n 7 d 5n - 1 T 2 a + (Ci + C 2 )nTd n - 1 ). 

Proof. At each h-th step of algorithm, if there are more 
than one roots of the corresponding system in the positive 
orthant (the cost of estimating this is C2, we compute the 
corresponding partial quotients lh = {lh,i, ••> lh,n)> where 
C- (hh,i) < o~h (the cost of estimating this is Ci Then, for each 
polynomial of the system, /, we perform the shift operation 
f(x\ + h, . . . ,x n + l n ), and then we split to 2™ subdomains. 
Let us estimate the cost of the last two operations. 

A shift operation on a polynomial of degree < d, by a 
number of bitsize cr, increases the bitsize of the polynomial 
by an additive factor nda. At the h step of the algorithm, 
the polynomials of the corresponding system are of bitsize 
0(t + ndy^ j _ 1 (Th), and we need to perform a shift opera- 
tion to all the variables, with number of bitsize an+i- The 
cost of this operation is Os(nd n r + n 2 d n+1 "J2h=i fffc )i an( ^ 
since we have n polynomials the costs becomes OB{n 2 d n r + 
n 3 d n+1 Ysk=i <Jfe )' The resulting polynomial has bitsize 0(t+ 

nd J2ktl °k). 

To compute the cost of splitting the domain, we proceed 
as follows. The cost is bounded by the cost of perform- 
ing n2 n operations f{xi + 1, . . . , x n + 1), which in turn is 
B (nd n T + n 2 d n+1 Yllt\ CT fe + n 2 d n+1 ). So the total cost 
becomes & B {2 n n 2 d n T + 2 n n z d n+1 Yllt\ <^k)- 

It remains to bound X]fc=i cr fc- If o" is a bound on the 
bitsize of all the partial quotients that appear during and 
execution of the algorithm, then J^fcii ° k = O(ho). 



Moreover, h < #(T) = 0(n 2 rd 2n - x ) (lfim.RTI). and so 
the cost of each step is Os(2 n n 5 d 3 "ra). 

Finally, multiplying by the number of steps (lem. 15. 1[) we 
get a bound of d B (2 n n 7 d 5n - 1 r 2 a). 

To derive the total complexity we have to take into ac- 
count that at each step we compute some partial quotients 
and and we count the number of real root of the system 
in the positive orthant. Hence the total complexity of the 
algorithm is d B (2 n n 7 d Sn ~ 1 r 2 a + (& + C 2 )nTd n ~ 1 ). □ 

In the univariate case (n — 1), if we assume that © holds 
for real algebraic numbers, then the cost of Ci and C 2 is 
dominated by that of the other steps, that is the splitting 
operations, and the (average) complexity becomes B (d 3 r) 
and matches the one derived in [18] (without scaling). 

5.3 Further improvements 

We can reduce the number of steps that the algorithm 
performs, and thus improve the total complexity bound of 
the algorithm, using the same trick as in |18| . The main 
idea is that the continued fraction expansion of a real root 
of a polynomial does not depend on the initial computed 
interval that contains all the roots. Thus, we spread away 
the roots by scaling the variables of the polynomials of the 
system by a carefully chosen value. 

If we apply the map (sri, ■ ■ ■ , x n ) 1— > (x\/2 l , . . . , x n /2 l ), to 
the initial polynomials of the system, then the real roots are 
multiply by 2 l , and thus their distance increase. The key 
observation is that the continued fraction expansion of the 
real roots does not depend on their integer part. Let ( be 
the roots of the system, and 7, be the roots after the scaling. 
It holds 7 = 2 l C,. From [8] it holds that 

- log Y[ MO < 2nrd 2n - 1 \g{nd 2n ) + 2nd" \g(nd 2n ), 
and thus 

-i og nA I (7)=-iog2 M n mo 

< (2nrd 2n - 1 + 2nd n ) \g{nd 2n ) - RL 

If we choose £ = 2nd n ~ 1 (d + r)lg(nd n ) and assume that 
R — d n which is the worst case, then — log Ilcev ^(7) = 0. 
Thus, following the proof of Lem. 15.11 the number of steps 
that the algorithm is 0(nd n ). 

The bitsize of the scaled polynomials becomes 0(n 2 d" +1 + 
n 2 d n r). The total complexity of algorithm is now 

d B (2 n n 5 d 3n+1 a + 2 n n 5 d 3n r + nd n (d+C 2 )), 

where a the maximum bitsize of the partial quotient appear 
during the execution of the algorithm. If we assume that 
(O holds for real algebraic numbers, then a — 0(1). Notice 
that in this case, when n = 1, the bound becomes B (d 3 r), 
which agrees with the one proved in [18| . 

The discussion above combined with Prop. I5.2l lead us to: 

Theorem 5.3. The total complexity of the algorithm is 
B (2 n n 5 d 3n+1 a + 2 n n 5 d 3n T + nd n (d+C 2 )). 

6. IMPLEMENTATION AND EXAMPLES 

We have implemented the algorithm in the C++ library 
realroot of MathemagixQi which is an open source effort 

2 http : //www . mathemagix . org/ 
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Figure 3: Isolating boxes of the real roots of (£1). 

that provides fundamental algebraic operations such as al- 
gebraic number manipulation tools, different types of uni- 
variate and multivariate polynomial root solvers, resultant 
and GCD computations, etc. 

The polynomials are internally represented as a vector of 
coefficients along with some additional data, such as a vari- 
able dictionary and the degree of the polynomial in every 
variable. This allows us to map the tensor of coefficients 
to the one-dimensional memory. The univariate solver that 
is used is the continued fraction solver; this is essentially 
the same algorithm with a different inclusion criterion, the 
Descartes rule. The same data structures is used to store the 
univariate polynomials, and the same shift/contraction rou- 
tines. The univariate solver outputs the roots in increasing 
order, as a result of a breadth-first traverse of the subdivi- 
sion tree. In fact, we only compute an isolation box for the 
smallest positive root of univariate polynomials and stop the 
solver as soon as the first root is found. Our code is tem- 
plated and is efficiently used with GMP arithmetic, since 
long integers appear as the box size decreases. 

The following four examples demonstrate the output of 
our implementation, which we visualize using AxeiQ. 

First, we consider the system /1 = fa = (Si), where 
/1 = x 2 + y 2 — xy — 1, and f 2 = Wxy — 4. We are looking 
for the real solutions in the domain I = [—2,3] x [—2,2], 
which is mapped to R+, by an initial transformation. The 
isolating boxes of the real roots can be seen in Fig. [3] 

In systems (£2), (S3), We multiply /1 and f 2 by quadratic 
components, hence we obtain 

, J h = + 2x 2 y 2 - 2x 2 + y 4 — 2y 2 - x 3 y - xy 3 + xy + 1 
1 21 \ f 2 = 20x 3 y ~ 10x 2 y 2 - Wxy 3 ~ 8x 2 + 4xy + 4y 2 

and 

{/1 = 10x 2 y - Wxy 3 - 4x + Ay 2 
f 2 =x 4 - 2x 2 y - 2x 2 + y 2 x 2 - 2y 3 - y 2 - x 3 y+ 
+2xy 2 + xy + 2y+l 

The isolating boxes of this system could be seen in Fig. [4] 
Notice, that size of the isolation boxes that are returned in 
this case is considerably smaller. 

Consider the system (£4), which consists of /1 = x 4 — 
2x 2 — y 4 + 1 and f 2 , which is a polynomial of bidegree (8,8). 
The output of the algorithm, that is the isolating boxes of 
the real roots can be seen in Fig. [4] 

3 http : //axel . inria.fr 



Figure 4: Isolating boxes of the real roots of the system Left: (E2), Middle: (S3), Right: (£4). 



System 


Domain 


Iters. 


Subdivs. 


Sols. 


Excluded 


Ei 


[0, 10f 


53 


26 


4 


25 


S 2 


[-2,3]" 


263 


131 


12 


126 


S3 


[-2,3]" 


335 


167 


8 


160 


s 4 


[-3,3]" 


1097 


548 


16 


533 



Table 1: Execution data for Si, S2, S3, S4. 

One important observation is the fact the isolating boxes 
are not squares, which verifies the adaptive nature of the 
proposed algorithm. 

We provide execution details on these experiments in Ta- 
bled Several optimizations can be applied to our code, but 
the results already indicate that our approach competes well 
with the Bernstein case. 
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